Decoupling the coupled DGLAP evolution equations: an analytic solution to pQCD 



Martin M. Block 

Department of Physics and Astronomy, Northwestern University, Evanston, IL 60208 

Loyal Durand 

Department of Physics, University of Wisconsin, Madison, WI 53706 

Phuoc Ha 

Department of Physics, Astronomy and Geosciences, Towson University , Towson, MD 21252 

Douglas W. McKay 
Department of Physics and Astronomy, University of Kansas, Lawrence, KS 66045 

(Dated: April 12, 2010) 

Using Laplace transform techniques, along with newly-developed accurate numerical inverse 
Laplace transform algorithms [l], U, we decouple the solutions for the singlet structure function 
F B (x,Q 2 ) and G(x,Q 2 ) of the two leading-order coupled singlet DGLAP equations, allowing us to 
write fully decoupled solutions: 

F s (x,Q 2 ) = F s (F s0 (x),G (x)), 
G(x,Q 2 ) = G(F a0 (x),G (x)). 

Here J-" s and Q are known functions — found using the DGLAP splitting functions — of the functions 
Fsq(x) = F s (x,Qo) and Gq(x) = G(x,Qq), the chosen starting functions at the virtuality Qq. As 
a proof of method, we compare our numerical results from the above equations with the published 
MSTW LO gluon and singlet F a distributions [jj, starting from their initial values at Q 2 , = 1 GeV 2 . 
Our method completely decouples the two LO distributions, at the same time guaranteeing that 
both distributions satisfy the singlet coupled DGLAP equations. It furnishes us with a new tool 
for readily obtaining the effects of the starting functions (independently) on the gluon and singlet 
structure functions, as functions of both Q 2 and Q 2 ,. In addition, it can also be used for non-singlet 
distributions, thus allowing one to solve analytically for individual quark and gluon distributions 
values at a given x and Q 2 , with typical numerical accuracies of about 1 part in 10 5 , rather than 
having to evolve numerically coupled integral-differential equations on a two-dimensional grid in 
x, Q 2 , as is currently done. 

Accurate knowledge of gluon distribution functions at small Bjorken x and large virtuality Q 2 plays a vital role in 
estimating QCD backgrounds and in calculating gluon- initiated processes, and thus in our ability to search for new 
physics at the Large Hadron Collider. 

The gluon and quark distribution functions have traditionally been determined simultaneously by fitting experi- 
mental data on neutral- and charged-current deep inelastic scattering processes and some jet data over a large domain 
of values of x and Q 2 . The distributions at small x and large Q 2 are determined mainly by the proton structure 
function i^ p (a;,Q 2 ) measured in deep inelastic ep (or j*p) scattering. The fitting process starts with an initial Qq, 
typically less than the square of the c quark mass, m 2 « 2 GeV 2 , and individual quark and gluon trial distribu- 
tions parameterized with pre-determined shapes, given as functions of x for the chosen Qq. The distributions are 
then evolved numerically on a two-dimensional grid in x and Q 2 to larger Q 2 using the coupled integral-differential 
DGLAP equations typically in leading order (LO) and next-to- leading order (NLO), and the results used 

to predict the measured quantities. The final distributions arc then determined by adjusting the input parameters 
to obtain a best fit to the data. This procedure is very indirect in the case of the gluon: the gluon distribution 
G(x,Q 2 ) = xg(x,Q 2 ) does not contribute directly to the accurately determined structure function F^ix, Q 2 ), and is 
determined only through the quark distributions in conjunction with the evolution equations, or at large x, from jet 
data. For recent determinations of the gluon and quark distributions, see 0, 0-Q3 ■ 

In the following, we will summarize our analytic method for determining the singlet structure functions F s (x,Q 2 ) 
and G(x,Q 2 ) directly and individually, using as input F s q(x) = F s (x,Qq) and Gq(x) = G(x,Qq), where Qq is 
arbitrary, with the guarantee that each individually satisfies the coupled DGLAP equations. The method can be 
extended simply to embrace non-singlet functions, so that it can be used to find individual quark distributions. 
However, we will not pursue that goal in this communication. Instead, we give a numerical demonstration by using 
the LO MSTW Q F s0 and Go at Qq — 1 GeV 2 to generate their singlet structure functions and gluon distributions at 
large Q 2 . Because our basic solutions are analytic, we readily obtain numerical accuracies for F s (x, Q 2 ) and G(x, Q 2 ) 
of better than 1 part in 10 5 for all Bjorken-x and virtuality Q 2 considered. 
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Our approach uses a somewhat unusual application of Laplace transforms [TT], [T2| , in which we first introduce the 
variable v = ln(l/a;) into the coupled DGLAP equations, then Laplace transform these coupled integral-differential 
equations in v space to obtain coupled homogeneous first-order differential equations in the variable Q 2 . The 
parameters of these transformed equations are known functions of s, the Laplace-space variable. These equations are 
then solved analytically. Finally, using fast and accurate numerical inverse Laplace transform algorithms [l], Q , we 
transform the solutions back into v space, and, finally, into Bjorken x-space, i.e., F s (x,Q 2 ) = J 7 s (F s o(x),Go(x)) and 
G(x,Q 2 ) = G(F s q(x),Gq(x)), where the functions T and Q are determined by the splitting functions in the DGLAP 
equations. 

Our method can be generalized to NLO, but for brevity, we will limit ourselves to LO in this paper. We write the 
coupled LO DGLAP equations [HI El as 

^ dF s ,~2n _ A1 p f r#\ , 16 „, n2u ,l^,16„ f 1 (F s {z,Q 2 ) F s (x,Q 2 )\ dz 



~r f Fs ^ Q2) i 1 + f ) $ + 2nfX f G(z ' Q2) i 1 - 2 - z + 2 5) $' (1) 

47r 9G ^,g 2 ).^^ G (,,g 2) + i2 G (,,Q 2 )inl^: + i2^V^Q 2 ) Gfcg^ & 



s (Q 2 )dlnQ 2 



X 



12x1 G(z,Q')[--2 + ---)— + - I F s {z,Q 2 ) [ 1+ (1- -J J — . (2) 



2\ dz 



Here a s (Q ) is the running strong coupling constant, given in LO by 

47T 

(11 - fn/) In(Q 2 /A2 / ) 

where A„, is fixed so that A 5 reproduces a s (M|), and the other A's (A4 and A3) arc adjusted so that a s is continuous 
across the boundaries Q 2 = M 2 and M 2 , respectively, where M& and M c are the masses of the b and c quarks. 

We now examine the last two terms of line 1 in Eq. (TT]) and rewrite them, introducing the variable changes v = 
ln(l/.x), w = ln(l/z), and the notation F s (v,Q 2 ) = F s ( e - V ,Q 2 ), G(v,Q 2 ) = G(e~ v ,Q 2 ) as 

Y F 8 (v, Q 2 ) He" - 1) + H jf fa™' ~ F ^ v ' Q 2 ) eV ~ W ) eV -l _ I dw 

= y£^KQ 2 )ln(l-e-(«- >) dw. (4) 

where the final result — the last line in Eq. Q — is found by replacing the upper limit v in integral of line 1 of Eq. (TJ| 
by v — e, carrying out the integrals, doing a partial integration and then taking the limit as e — > 0. Similarly, we find 
for the last two terms of line 1 in Eq. , that 



12G(v, Q 2 ) ln(e» - 1) + 12 J (d(w, Q 2 ) - G(v, Q 2 )e v ~ w ) — _ { dw 

-(ly.Q^lnfl-e-^-^) dw. (5) 
1 V / 



" dG { 
dw 
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We next rewrite Eq. ([T]) and Eq. ([2]), in terms of the new variable v, as 



4tt dF a 
a s (Q 2 )d In Q' 



■(v,Q 2 ) = AF s (v,Q 2 ) + 



16 r dF s 
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3 Jo 



F, 



3 Jg dw 
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«,Q 2 ) + 12 y ^( w ,g 2 )ln(l- e -( u -^J <2w 



(6) 
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+ 12 y G(w, Q 2 ) (l - 2e-^ v ~ w) + e' 2 ^-^ - e^"™)) dw 
+ ^ jT F s (w, Q 2 ) (l + (l - e-(— )) dw. 



(7) 



All of the integrals in Eq. © and Eq. ([7]) are convolutions. Introducing Laplace transforms allows us to factor these 
integrals, since the Laplace transform of a convolution is the product of the Laplace transform of the factors, i.e., 



= C[F[v];s}xC[H[v];s}. 



r rv 

C / F[w]H [v — w] dw; s 
Jo 

Defining the Laplace transforms 

f( s ,Q 2 ) ee £[f s (v,Q 2 );s] , g(s,Q 2 )=£[G(v,Q 2 );s] 

and noting that 
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(8) 
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we can factor the Laplace transforms of Eq. © and Eq. Q into two coupled ordinary first order differential equations 
in Laplace space s with <5 2 -dependent coefficients. These can be written as 



dlnQ 2 

dg 



(s,Q 2 ) 

(s,g 2 ) 



<91nQ 2 

whose coefficients $ and are given by 
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(11) 
(12) 

(13) 
(14) 
(15) 
(16) 



where ip{x) is the digamma function and — 0.5772156 ... is Eulcr's constant. 

The solution of the coupled equations in Eq. (fTTj) and Eq. (jl"2j) in terms of initial values of the functions / and 
g, specified as functions of s at virtuality Qq, is straightforward. The Q 2 dependence of the solutions is expressed 
entirely through the function 



Q 2 

t(Q 2 ,Q 2 o) = ^- f a s (Q' 2 )d\nQ' 

47r Jc& 



(17) 
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With the initial conditions fo(s) = f(s, Qq) and go(s) = g(s, Q%), the solutions are 



/( S)T ) = fc//(s,r)/o(s) + kfg(s,T)g (s), 
g(s,r) = k gg (s,T)g (s) + k gf (s,T)f (s), 



where the coefficient functions in the solution are 



k ff{s,r) 

kfg(S,T) 

kgg( s , T ) 

kg f (a,T) 



= pi(*/W+* 9 W) 



cosh ( -R(s] 



sinh (§i?(s)) 
W) 
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(18) 
(19) 



(20) 
(21) 
(22) 
(23) 



with i?( S ) = t/ (*/(*) - *»(*)) + 46/(5)63(5). 



Let us now define four kernels K FF , K FG , K GF and K GG , the inverse Laplace transforms of the k's, i.e., 

K ff (v,t) = £~ 1 [fc//(s,r);u], K fg (v,t) = £~ 1 [k fg (s 1 T);v], 
K gg (v,t) = C^ 1 [k gg {s,T)- 1 v\, K gf (v,t) = C^ 1 [k gf {s,T);v\. 



(24) 
(25) 



It is evident from Eqs. (fTT)) . ([2"T]) . and (j2"3")) that IffG and K GF vanish for Q 2 = Qq where t(Q 2 ,Qq) = 0. It can also 
be shown without difficulty that for r = 0, K FF (v, 0) = Kqq{v, 0) = S(v). 

The initial boundary conditions at Qq are given by F s q(x) = F s (x, Qq) and Gq(x) = G(x, Qq)- In v-space, 
F s q(v) = F s0 (e~ v ) and Go(t>) = Go(e~") are the inverse Laplace transforms of fo(s) and go(s), respectively, i.e., 



F s0 (v) = C- 1 [fo(s);v}andGo(v) = C- 1 [g (s);v]. 
Finally, we can write our decoupled solutions in w-space in terms of the convolution integrals 

F s (v,Q 2 ) = 

G(v,Q 2 ) = I K GG (v-w,T(Q 2 ,Q 2 ))G (w)dw + 



K FF (v-w,T(Q 2 ,Q 2 Q))F s0 {w)dw+ / K FG (v - w,t(Q z ,Q 2 ))G (w) dw 



(v-w,t(Q',Q 2 ))F s0 (w) dw. 



(26) 

(27) 
(28) 



Noting again that v = ln(l/x), in the usual variables — Bjorken-x and virtuality Q 2 — we readily find the desired 
decoupled F s (x,Q 2 ) and G(x,Q 2 ) from the above decoupled solutions for F s (v,Q 2 ) and G(v,Q 2 ), requiring only a 
knowledge of the initial values F s0 (x) and Go(x) at Qq. 

For non-singlet distributions F ns (x,Q 2 ), such as the difference between the u and d quark distributions, 
x [u(x, Q 2 ) — d(x, Q 2 )] , we can schematically write the logarithmic derivative of F ns as the convolution of F ns (x, Q 2 ) 
with the non-singlet splitting function lC ns (x) (using the convolution symbol ®), i.e., 



4^ dF n 



-(x,Q 2 ) = F„ s ®K f 



(29) 



a s (Q 2 )d\n(Q 2 ) 

After changing to the variable v = ln(l/x) and going to Laplace space s, we find the simple solution 

f„ s (s,T) = e T *" s(s) /n,o(s), where $„ s (s) = £ [ e -"/C ni (u); s] and t ns (v) = K na (e"") . (30) 

Thus we can find any non-singlet solution in w-space, using the non-singlet kernel K ns (v) = C^ 1 [e 1 "*" 3 ^; u] , by 
employing the Laplace convolution relation 



F ns (v,Q 2 



K ns (v - w,T(Q 2 ,Ql))F ns0 (w)dw. 



(31) 
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For brevity, we will not pursue the case of the non-singlet solution any further here except to note that in LO, the 
&ns( s ) m Eq. (f5U|) is identical to ^/(s) defined in Eq. (|13p . but will concentrate instead on the more difficult case of 
F s and G. 

As an example of the application of this technique, we will compare the x-space singlet distribution function 
F.(x, Q 2 ) calculated from Eq. {27]) starting from the MSTW initial conditions at Q\ = 1 GeV 2 , with the LO MSTW 
[3| distributions. We also will compare separately their G(x,Q 2 ) with the results found from Eq. (|28]l. We evaluate 
the kernels Kff(u,t), Kfg( u t1~) m Eq. (|27|) and Kgg{i>' 1 t) 1 Kqf(u,t) in Eq. (|28|) numerically, using powerful new 
inverse Laplace transformation algorithms p], ■ In order to insure continuity across the boundaries Q 2 = M 2 and 
M 2 , we will first evolve from Qq = 1 GeV 2 (the MSTW Q\ value) to M 2 and use our evolved values of G(x, M 2 ) and 
F s (x, M 2 ) for new starting values G (x) and F s0 (x). We will then evolve to M 2 , repeating the process, thus insuring 
continuity of G and F s at the boundaries where rif changes. We use the MSTW values M c = 1.40 GeV, M c = 4.75 
GeV, a s (l GeV 2 ) = 0.6818 and a s {M 2 z ) = 0.13939. 

In Fig. [T] we show the results — in x-space — for LO G(x,Q 2 ) = xg(x,Q 2 ), for 4 values of Q 2 . The curves are 
the published MSTWjl LO gluon distributions, for Q 2 = 5, 20, 100 and M 2 GeV 2 , bottom to top. Since the 
MSTW collaboration |3[ started their evolution at Qq = 1 GeV 2 , we used F s o and Go constructed from their values 
at Qq = 1 GeV 2 in Eq. (f2"5)l . The dots are our results for LO G(x, Q 2 ) from Eq. (|2"5)l (converted to x-space), using the 
LO MSTW values for F s0 (x) and Gq(x). The agreement, over this large span of Q 2 , is quite striking. Our numerical 
accuracy has been investigated and is typically better than about 1 part in 10 5 at small Bjorkcn-x. The most serious 
disagreements between our calculated G and the MSTW curves are, at x= 10 -5 , 1.8% and 1.6% at Q 2 = 100 GeV 2 
and M§, respectively, which is approximately within their stated numerical errors. 

In Fig. [2] we show the results — in x-space — for the LO singlet distribution F s (x, Q 2 ), for 4 values of Q 2 , where again 
F s o and Go are the MSTW [|| values at their initial evolution value of Q 2 = 1 GeV 2 . The curves are the published 
MSTW LO singlet distributions, for Q 2 = 5, 20, 100 and M 2 GeV 2 , bottom to top. The dots are our results for 
LO F s (x,Q 2 ) from Eq. (|2"7) (converted to cc-space), using the LO MSTW values for F s0 (x) and Go(x). Again, the 
agreement is excellent over the entire range of Bjorken-x and virtuality Q 2 . The most serious disagreements between 
our calculated F s and the MSTW curves are, at x= 10~ 5 , 2.0% and 1.7%, at Q 2 = 100 GeV 2 and M§, respectively. It 
is clear from Fig. Q]and Fig. [5] for G s and F s — over the enormous virtuality and x span — that our analytic solutions 
of Eq. (f2"T|) and Eq. (|2"5)) arc correct. The numerical values were evaluated using Mathematica (l3| . 

In conclusion, we have constructed decoupled analytical evolution equations for F s (x,Q 2 ) and G(x,Q 2 ) from the 
coupled LO DGLAP equations. These require only a knowledge F s q(x) and Gq(x), the initial values of F s and G at 
the starting value Q\ for the evolution, to calculate F s (x,Q 2 ) and G(x,Q 2 ). The same procedures can be used for 
non-singlet distributions, allowing one to obtain analytic solutions for individual quark distributions, as well as for the 
gluon distribution, avoiding the necessity for numerical solutions of the coupled DGLAP on a two-dimensional grid 
in x, Q 2 . In essence, in a program such as Mathematica, we could define a function for each quark and gluon, and by 
inputting the desired x and Q 2 , simply evaluate it, accurately and rapidly. 

We demonstrated numerically that the method gives agreement with published MSTW Q LO values of G{x,Q 2 ) 
and F s (x, Q 2 ) over an enormous range of x and Q 2 . The accuracy obtained using our analytic solution and fast new 
algorithms for performing inverse Laplace transforms 0, [|J was better than 1 part in 10 5 . In the future, as well as 
evaluating non-singlet distributions and the NLO singlet case, we will evaluate F s q(x) and Go(.t) in both LO and 
NLO, from a fit to small x experimental data for the structure function F^ix^Q 2 ), in order to analytically obtain 
accurate values of G(x, Q 2 ) that are directly tied to experiment, which are needed for the LHC. 
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work were done. P. Ha would like to thank Towson University Fisher College of Science and Mathematics for travel 
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FIG. 1: The LO MSTW gluon distribution, G(x,Q 2 ) = xg(x,Q 2 ), for Q 2 = 5 ,20, 100 and M% GeV 2 . The published 
MSTW 0| curves are for Q 2 = 5, 20, 100 and M 2 GeV 2 , bottom to top. The dots are the evolution results for LO G(x,Q 2 ) 
from Eq. ([28} (converted to a;-space), using the LO MSTW values for F s o(x) and Go (as), where Qq = 1 GeV 2 . 
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FIG. 2: The LO MSTW singlet distribution, F s {x,Q 2 ), for Q 2 = 5 ,20, 100 and M§ GeV 2 . The published MSTW 
curves are for Q 2 = 5, 20, 100 and M 2 GeV 2 , bottom to top. The dots are the evolution results for LO F s (x, Q 2 ) from Eq. (|27|) 
(converted to x-space), using the LO MSTW values for F s0 (x) and Go (a;), where Ql = 1 GeV 2 . 



